Hopping dynamics of interacting polarons 
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We derive an effective cluster model to address the transport properties of mutually interacting 
small polarons. We propose a decoupling scheme where the hopping dynamics of any given particle is 
determined by separating out explicitly the degrees of freedom of its environment, which are treated 
at a statistical level. The general cavity method developed here shows that the long-range Coulomb 
repulsion between the carriers leads to a net increase of the thermal activation barrier for electrical 
transport, and hence to a sizable reduction of the carrier mobility. A mean-field calculation of 
this effect is provided, based on the known correlation functions of the interacting liquid in two and 
three dimensions. The present theory gives a natural explanation of recent experiments performed in 
organic field-effect transistors with highly polarizable gate dielectrics, and might well find application 
in other classes of polaronic systems such as doped transition-metal oxides. 



I. INTRODUCTION 

In recent years, the development of organic electron- 
ics has triggered a strong effort towards the understand- 
ing of charge transport in organic field-effect transistors 
(OFETs) i In such devices, the carriers induced by a gate 
potential move at the interface between an organic semi- 
conductor and a dielectric. Unlike their inorganic coun- 
terparts, such as Si MOSFETs, the transport properties 
in OFETs are dominated by the weak transfer integrals 
between the molecular constituents of the organic ma- 
terial: the Van der Waals inter-molecular bonding leads 
even in pure crystalline samples to extremely narrow elec- 
tronic bands, making such systems very sensitive to in- 
teractions. 

Recently, a systematic study of rubrene-based single- 
crystalline OFETs fabricated using gate materials of in- 
creasing dielectric polarizability has revealed that, in the 
case of high-K dielectrics, the dominant limiting mech- 
anism of electron transport originates from the cou- 
pling with the polar phonons at the organic-dielectric 
interface.— This phenomenon, which in wide-band inor- 
ganic semiconductors 3 -* 4 ^^ and in graphene^ only leads 
to minor modifications of the electron mobility, can be 
so effective in organic semiconductors that it leads to 
polaronic self-localization of the carriers on the scale of 
one or few molecules. As a consequence, the mobility 
is strongly suppressed and becomes thermally activated, 
being due to the incoherent hopping of small polarons on 
the molecular lattice. 

Because of the increased capacitance of the devices, 
the use of high-K dielectrics also has a second interesting 
consequence, as it allows the injection of sufficiently large 
charge densities, such that the electrons can no longer be 
considered as non-interacting carriers^ Indeed, concen- 
trations of the order of 0.1 carriers/molecule and above 
have been reached in rubrene devices using Ta 2 05 as a 
gate material (dielectric constant e s = 25). The current- 



voltage characteristics of such devices exhibit strong de- 
viations from linearity that cannot be explained in terms 
of independent carriers, and have been ascribed to the 
onset of electron-electron interactions^ 

The aim of this work is to establish a theory for the 
density-dependent transport properties of mutually in- 
teracting small polarons in the hopping regime. Al- 
though the present derivation is motivated by the physics 
of organic-dielectric interfaces, the problem itself is suf- 
ficiently general to find application in other polaronic 
systems such as transition-metal oxides 1 ^ and possi- 
bly oxide-oxide interfaces 1 ^ and organic charge transfer 
interfaces^ As will become clear in the following, how- 
ever, OFETs are ideal systems for the observation of the 
many-body effects studied here, for two reasons. First, 
the effect of interactions on the polaronic hopping rates 
can be comparatively large in organic semiconductors, 
where small polarons can exist with activation energies 
that are generally smaller than in oxides. Secondly, and 
most importantly, in such devices the carrier concentra- 
tion can be varied accurately by tuning the gate voltage 
without the need of chemical substitutional^ 3 - thus pro- 
viding a reliable and unambiguous procedure to disentan- 
gle many-body effects from the properties of individual 
polarons. 

The present paper provides a detailed derivation of the 
theory that was used in Ref]^ to explain the current- 
voltage characteristics of rubrene-based high-K OFETs, 
as well as its generalization to three-dimensional systems. 
We start by introducing an effective cluster model that al- 
lows to separate the dynamics of any given electron from 
the remaining carriers in the system, which are treated as 
an external statistical environment. The mobility is then 
obtained by assuming a succession of incoherent hopping 
events, from a suitable statistical average of the inter- 
molecular hopping rates. The theory is applied to the cal- 
culation of the density-dependent transport properties of 
small polarons interacting through long-range Coulomb 



2 



potentials. It is shown that the mutual interactions be- 
tween carriers gives rise to a net increase of the polaronic 
thermal activation barrier, and consequently to a sizable 
reduction of the mobility. A mean-field calculation is 
performed, providing a closed-form analytical expression 
for the mobility based on the known pair distribution 
function of the interacting liquid. 

The paper is organized as follows. In Sec. |n]we carry 
out explicitly the separation between the cluster degrees 
of freedom and those of the environment. In Sec. IIIII 
we focus specifically on a two-site cluster, which is ap- 
propriate in the small polaron limit. We derive a gen- 
eral formula for the hopping mobility of interacting small 
polarons and provide a simplified expression describing 
carrier-carrier correlations in the mean-field approxima- 
tion. In Sec. IIVI the theory is applied to the problem of 
polarons interacting via the long-range Coulomb repul- 
sion. The results are discussed in relation to the trans- 
port properties of OFETs and of other known polaronic 
systems in Sec. [V] 



II. EFFECTIVE CLUSTER MODEL 

Let us consider a system of electrons moving on a lat- 
tice (defined by the lattice vectors Rj) which interact 
with each other and with the lattice vibrations, as de- 
scribed by the following Hamiltonian: 

H = H t + H p h + Hi n t + H e ^ p h- (1) 
The first term 

H t = -tJ2 4*3 (2) 

<ij> 

is the tight-binding Hamiltonian for free electrons, with 
t the hopping integral and Cj, c\ the corresponding de- 
struction and creation operators for electrons. The sec- 
ond term 

describes local (dispersionless) phonons of frequency 
luq = y/k/M, where Xj is the local lattice displacement, 
k the spring constant and M the mass. The term 

H mt = ]^^n l Vi n J (4) 

i,j 

is the electron-electron (e-e) interaction where rii = cjcj 
is the electron density operator at site i, and Vij = 
V(Ri — Rj) a generic density-density interaction poten- 
tial. Finally, 

H e -ph = ^2 riigijXj (5) 



is the electron-phonon (e-ph) interaction where gij = 
g(Ri — Rj) is a non-local density-displacement coupling. 
We shall not give here the precise form of and Vij , 
which is not needed for the general formalism developed 
in the following paragraphs. This will be provided later, 
when dealing with specific examples. 

A. Derivation of the cluster Hamiltonian 

We now derive an effective cluster model neglecting 
the second term in Eq. ((3]), which is valid at tempera- 
tures T > luq. The phonon kinetic energy term will be 
reintroduced in Section II B to treat the polaron hopping 
dynamics in the semiclassical adiabatic approximation. 

Our starting point to evaluate the hopping mobility of 
the system described by Eq. (JTJ is to divide the lattice 
into a cluster (c) of finite size, in which we allow for 
quantum electron hopping, and a remaining part (c) in 
which the dynamics of the electrons is neglected. This 
separation is enforced by keeping a finite transfer integral 
t only for electrons within the cluster (c), while setting 
t = in (c). It is then natural to rewrite the model Eq. 
([I]) by singling out the terms which explicitly contain 
electronic variables in (c), that we denote as H,: 

H.=H^+H% + Hi% h + (6) 

so that 

H = H ph + H% t +Hf} ph +H.. (7) 

In Eqs. © and ([7]) the labels (c) and (c) indicate that 
the sums over electronic variables are restricted respec- 
tively to the cluster or the environment, and (c, c) stands 
for interactions among electrons belonging to the two dif- 
ferent sub-systems. It can be noted that all terms in Eq. 
(J7J) commute with each other. 

The proposed separation scheme is formally equivalent 
to the one used by Pardee and Maha n 14 ' 15 to describe 
electrical conduction in solid electrolytes. The justifica- 
tion in that case follows from the large masses of the 
charge carriers, that are mobile ions. As a result, the 
collective rearrangement of the particles in reaction to a 
given hopping event is much slower than the hopping pro- 
cess itself, so that the positions of the carriers in the en- 
vironment can effectively be regarded as static variables 
during the local dynamical evolution. In our case this 
decoupling is justified due to the exponential suppres- 
sion of the carriers' hopping rate associated to polaronic 
self-localization. 

To derive an effective cluster model it is useful to in- 
troduce the following reduced density matrix: 

Pr = |tr ( ^ ) tr (E) e- H /' £ - T n( c ) ( 5(y i - ^ffyXj) (8) 

3 

where n^ c ^ indicates the product over the cluster elec- 
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irons, the trace symbols are defined as 



tr 



(c) 



(..) = E (■■) 



rii ,i£(c) 

tr(p/i)(..) = / UidXi(..), 



(9) 



(10) 



and Z = tr( c )tr( p/l )tr( g )e~ ff / fc£iT . In Eq. (|9|) we trace 
over all the electronic degrees of freedom which do not 
belong to the cluster. The trace over phononic variables 
in Eq. (|10|) is performed by assigning the value of the 
phonon-induced external fields 



Y i = ^2g ij X j t6(c). 



(11) 



These are the phononic collective variables which act on 
each site of the cluster, through the non-local e-ph inter- 
action Eq. (0. The two steps described above are now 
explicitly carried out. 

a. Tracing out the phonons. The term H m in Eq. © 
depends on the phonons only through the variables {Y{\. 
We can therefore rewrite Eq. ([8]) as 



o - -tr, Ntr, -e-( H '« y *}) +i? £*)/ fesT 

Pr — ^^ L (c) lL (ph) t: 



(12) 



x 



nr^-£^)e-^+^ 



)/k B T 



where with H,({Yi}) we indicate the explicit dependence 
of this term on the collective phonon variables. The trace 
over the original phonons Xj in Eq. (|12p can be per- 
formed by introducing the integral representation of the 
S function 

8<y t - = / ^e**«-£,«,*i>. (13) 

3 

Performing the gaussian integrals over {X{} and over 
{wj} we obtain 



p r oc tr (g )e 



-(H.{Y)+H^ t +H^ f )/k B T 



(14) 



where H^ff represents the effective Hamiltonian resulting 
from the trace over phonons. It can be expressed as 



1 



(15) 



with 



1 1 (c) 

D *i = t[9%3 - rJ2^ 2 U9cWh (16) 



l.k 



H '% = tE^ 2 ]^+EE G ^ 



(<0 (e) 



(c) 



(18) 



In the above equations we have introduced the symbol 
[g~ 2 ] to denote the inverse of the matrix [g 2 ] in the cluster 
sub-space. As can be seen from Eq. (fl5]) . integrating out 
the phonon variables has led to an effective attraction 
Dij between the (c) electrons, whose form is given by 
Eq. (fT6|) . Similarly, equations (fi"7]) and (fl"8|) describe 
the effective interactions arising between the (c) electrons 
and the collective variables Yf. 

Adding the phonon-mediated interaction of Eq. (|15p to 

the bare electron-electron term H^ t in Eq. ([7]) yields the 
following screened interaction between the environment 
electrons 



H 



- E m (Vij 



D 4 



(19) 



1,3 



The reduced density matrix can be finally expressed as 

p r oc t m e-( H '^ + "^ +H '$fV k ° T . (20) 

Before moving on to the integration of the environment 
electrons it is useful to comment on the physical meaning 
of the two different contributions to the phonon-induced 
screening in Eq. (fTB]) . The first term, which leads to the 
effective potential 



V« = Vii 



[9% 
k 



(21) 



represents the ability of the polarizable medium to par- 
tially screen the electron-electron interaction. For ex- 
ample, starting from the bare Coulomb potential Vij = 
e 2 / 1 {^ooRij) and an electron-phonon interaction gij of the 
Frohlich type, it is shown in Appendix [Bl that the inclu- 
sion this term yields Vij — e 2 /(e s i?y), which correctly 
reproduces the static screening response of a bulk polar 
dielectric. 

The second term in Eq. (flBj) is a cavity field which 
arises due to the constraints in Eq. (fTTj) , because not 
all of the phonons have been integrated out. It can be 
viewed as the part of phonon screening that is missing 
due to the existence of the cluster. Since it involves the 
product of two matrices [g 2 ]u [g 2 ]kj , which decays faster 
than the direct screening [g 2 ]ij itself, this term becomes 
negligible when the cluster size is smaller than the av- 
erage interparticle distance. This cavity correction can 
therefore be neglected to lowest order in the electron con- 
centration, although its actual magnitude depends on the 
shape of the electron-phonon interaction g$j (for exam- 
ple, such cavity field is clearly absent in the limit of local 
e-ph interactions, i.e. g^ oc 5ij). On the other hand, if 
the cluster is enlarged to attain the size of the entire sys- 
tem, the two terms in Eq. (|16[) exactly cancel, and only 
the bare electron-electron interaction remains. 

b. Tracing out the electronic environment. The 
trace appearing in Eq. (|20p can be formally carried out 
by introducing two classical fields which couple linearly 
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to the cluster variables Yi and n», namely: 

(e) 

Vi = ^2 GijUj ; i e c (22) 

3 

(E) 

e,; = /^Vjjnj ; (23) 

3 

Such fields take into account the interactions between 
electrons in (c) and the cluster degrees of freedom, as con- 
tained explicitly in H'^ff and in the direct term H^ t . 
Substituting these definitions into Eqs. ^ and (TTT)) and 
regrouping terms in Eq. (|20p one obtains the following 
cluster Hamiltonian: 

0) x (c) (c) 

Hduster — ~t c]cj + o y^ n '^' n J + ^ n » e i 

<ij> i,j i 

(c) , (c) 

+ 1] % + rii) + - X 1 -j^'j- (24) 

Finally, by enforcing the definitions Eqs. (f22l) . (|23]) 
through the appropriate 6 functions, the reduced density 
matrix of the cluster can be expressed as a trace over the 
classical variables rji , ti of the environment 

Pr = J U^dmdete- 11 ^ /k B Tp ({£ . }) {7? . }) f (25) 
whose statistical distribution is 

(E) 

(s) 

x^-^V^)- (26) 

To summarize, Eqs. ([23]). ((25|) and (26]) describe a fi- 
nite cluster in which electrons mutually interact via the 
bare potential Vij , and are coupled to collective phonon 
variables Yi . The cluster degrees of freedom are also sub- 
ject to random fields rji and e; arising from the environ- 
ment electrons. Such fields are distributed, via Eq. (f2T))) . 
according to the equilibrium distribution of classical par- 

~ (c) 

tides interacting through the screened Hamiltonian H^ nt 
defined in Eq. HID. 

B. Ehrenfest dynamics of the cluster model 

Within the adiabatic regime, the carrier motion is 
constrained to follow the slow dynamics of the phonon 
coordinates . 10 ! 16 To determine the polaron mobility it is 
therefore necessary to treat explicitly the dynamics of 
the Xi that was neglected in the preceding Section. This 
can be done by introducing the semi-classical evolution 



of the lattice degrees of freedom through the following 
Ehrenfest equations 

MX, = -kXi - J2 9v ( n 3 (*)> ■ (27) 

3 

In the above equation the average of the electronic op- 
erators is taken at fixed Xi. To change to the cluster 
variables Yi we substitute Eq. (|2~7| into Eq. (jTTJ) for 
ie (c), leading to: 

(c) (E) 

M% = -kYi - J2l9ch <%•(*)) ~ I> 2 ]i>i>> (28) 

3 3 

where we have made explicit use of the assumption that 
the environment electrons do not evolve in time. Using 
Eq. (f2"2")) this can be rewritten as 

M% = -kY { - 5> c 2 k [{nj(t)) + rjj] . (29) 

3 

The collective phonon variables Yi are therefore subject 
to an external force which depends both on the instan- 
taneous electron density within the cluster and on the 
environment degrees of freedom through the fields rji. It 
is interesting to observe that in the present treatment, 
the frequency of the collective modes is equal to the bare 
phonon frequency ujq — y/k/M. The above Eq. (j2"9")l 
can equivalently be derived in a Hamiltonian formula- 
tion, by adding a kinetic term (2M) _1 J2if [9c]ij^i^j to 
Eq. with n, the momentum conjugate to Yi. 

ILL SMALL POLARON LIMIT 

The actual choice of the cluster size for practical cal- 
culations is dictated by the polaron properties, since 
it should be large enough to accomodate the electronic 
wavefunction involved in the hopping process. To keep 
the discussion simple and provide a physically signifi- 
cant example of the theory presented so far, we now 
focus specifically on the small polaron limit, where the 
electronic wavefunction collapses onto a single molecule. 
This situation is realized in systems with narrow elec- 
tronic bands, provided that the electron-phonon coupling 
is sufficiently strong. To be specific, this occurs when the 
energy of a polaron fully localized on a single molecular 
site, Ep = [g 2 ]n/2k, is larger than approximately half 
the free electron bandwidth, in which case a self-localized 
state becomes energetically more favorable than an ex- 
tended wave. The proper cluster in this case consists 
of two molecules — the initial (filled) site and the final 
(empty) site — and constitutes the basis for the theory 
of small-polaron transporti 10 ' 16 ' 17 We shall explicitly con- 
sider situations where the formation of bipolaronic states 
is ruled out by the presence of sufficiently strong repulsive 
interactions between the carriers ! 18 ' 19 Apart from this re- 
striction, the results obtained in this Section concerning 



5 



the effect of electron-electron interactions will be gener- 
ally valid regardless of the physical origin, and particular 
form, of gij , the only requirement being that the polarons 
are small. 



A. Two-site cluster 

It is shown in Appendix |A] that for a singly occupied 
two-site cluster the Hamiltonian Eq. (I24[) reduces to a 
spin-boson model, where the electronic degree of freedom 
plays the role of a pseudo-spin. Introducing the notation 
&z — n i — n 2 and cr x — c^c 2 + c^Ci we obtain 

H sb = -ta x + ^kQ 2 - -±={gQ + £)<t z . (30) 

The relative electronic occupation a z is coupled to the 
phonons through a single "interaction coordinate" 



Q 



V2 



9 



9_V2 
k 



V2 



(31) 



The first term in Eq. pip is the direct interaction 
with the collective phonons, and the second term orig- 
inates from the residual electron-phonon interaction of 
Eq. ([IS]), ryi and r] 2 being defined by Eq. ([22]). The 
parameter g is an effective electron-phonon coupling for 
the two-site cluster, defined through 



9 



[9 2 }n 



[g 2 



(32) 



The coupling with the environment electrons (c) also oc- 
curs via a single classical variable 



(£2 - ei) - y(?? 2 



(33) 



which takes into account the electronic repulsion on the 
two sites, corrected by the appropriate phonon mediated 
attractive terms Tfc. Using Eq. (|21|) , this can be rewritten 
as 



(5) 



(34) 



Such "local field" represents the energy difference be- 
tween the two-sites of the cluster in the presence of the 
potentials of the remaining electrons, screened by the lat- 
tice polarization. 



B. Adiabatic hopping 

In order to determine the polaron hopping rate, we 
now calculate the evolution of the dynamical variables 
a z and Q within the cluster in the presence of the lo- 
cal field £, which by assumption is fixed during the 
time of the hopping process. The electronic variable 



a z evolves quantum- mechanically through Eq. (|30|) . 
while the phonon collective variable Q is taken to evolve 
through the classical Ehrenfest equations 2 ^ Eq. (j2"9"j) . 
that reduce to 



M 



d 2 Q 
dt 2 



-kQ 



V2 



(35) 



where the average of the pseudo-spin is taken at a given 
configuration Q(t). 

As a further approximation, we estimate the electron 
transiti on pro bability within the adiabatic formulation 
of Refs. [l6lfl7l If the electron dynamics is faster than the 
motion of the phonons, the quantum variable a z is able 
to equilibrate at any given value of the classical Q. In 
this approximation the right-hand side of Eq. (|35p can 
be obtained from the derivative with respect to Q of the 
following adiabatic potential 



V ad (Q) = -k B Tlogtr a 



-H sb (Q,<?)/k B T 



(36) 



At sufficiently low temperature [lower than the barrier 
A(£) defined below], the adiabatic potential reads: 



V ad (Q) = \kQ 2 - V %v / 2 + .gQ) 2 /2 + t 2 



(37) 



In the polaronic regime, it has the double-well shape il- 
lustrated in Fig. [Th.. Within the adiabatic description, an 
electron at site 1 is associated to a phononic variable be- 
ing at the equilibrium point Q ~ —g/y/2k. For the occur- 
rence of a classical non-dissipative motion to the neigh- 
boring well, the initial kinetic energy of Q taken from a 
Maxwell distribution must exceed the relative maximum 
of V a d- This defines an energy barrier 

A(0 = Ap + | + gJ + O(t 2 ), (38) 

where Ap = ,g 2 /4/e — t is the activation barrier for inde- 
pendent adiabatic polarons,— The rate of electron hops 
per unit time from site 1 to site 2 then acquires a char- 
acteristic thermally activated behavior: 



w(0 



Wo 

2tt 



exp[-A(£)/fc B T] 



(39) 



We note that ui^ is the time it takes for the classical 
coordinate to pass from the initial point at Q ~ —gj\f2k 
to the final point at Q ~ g/V2k, and therefore corre- 
sponds to the natural timescale of the hopping process. 
As can be seen from Eq. (|39|) . in the hopping regime 
the time w(£) -1 between hopping events is exponentially 
longer than the hopping timescale uJq 1 . This fact al- 
lows to neglect the rearrangement of the environment 
electrons during a given hopping process, validating the 
cluster/environment separation scheme proposed in the 
present work. 

Finally, we remark that the adiabatic treatment devel- 
oped here is valid when the condition 



hu>o 

IT 



(g 2 /4k)k B T 



1/2 



(40) 



6 



^ ad 




'Q 




its own local field £. The calculation of the mobility 
therefore amounts to averaging the hopping rate over all 
the possible values of the local field £ through the ap- 
propriate distribution P(£). The mobility can then be 
written through Einstein's relation as 



(42) 



b) 
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FIG. 1: a) Adiabatic double-well potential V a d(Q) of Eq. (|37j) 
for a two-site cluster, in the absence (left) and in the presence 
(right) of electron-electron interactions. The electron position 
follows the dynamics of the phonon coordinate Q between the 
two minima of the potential. As pictorially represented in 
panel b), interactions between the carriers modify the shape 
of the double-well via the local field £, which measures the 
energy unbalance between the two sites due to the presence 
of the environment electrons. The grid represents the molec- 
ular lattice sites and the shaded area is the two-site cluster. 
The initial and final site of the carrier hop are indicated re- 
spectively by a black dot and a white circle. 



is meti 16 ' 17 In the opposite non-adiabatic regime, a 

hopping rate equivalent to the one given by Eqs. 

(f3"5| and ([55)) is obtained, although with a different 
prefactor-^2L22^3 



4TA, 



1/2 



cxp[-A(£)/fc B T] 



(41) 



and with Ap ~ (g 2 /4k). This has exactly the same de- 
pendence as Eq. ([39]) on the local field £ which embod- 
ies the effects of carrier-carrier interactions. It therefore 
appears that the theory developed here for interacting 
small polarons holds independently of the adiabatic/non- 
adiabatic character of the polaronic transport [Eq. ([40|) ]. 
provided that the appropriate prefactor is used in the 
hopping rate. 



C. Small polaron mobility 

To determine the mobility, we assume that the current 
flow occurs through a succession of incoherent hopping 
events. Each individual process is characterized by a rate 
of the form Eq. (f3~9~|) . which depends explicilty on the 
electronic environment of the hopping particle through 



where a is the length of the electron hop, which we take 
to be equal to the inter-molecular distance, and (w) is 
the statistical average 



(43) 



To find the statistical distribution that enters in Eq. 
([4"3"| we observe that if the system is sufficiently close 
to equilibrium, the value of the (static) local field £ is 
determined, via Eq. (|34p , by the positions of the en- 
vironment electrons prior to the hop. Correspondingly, 
P(£) follows, via Eq. (f26|) . from the equilibrium distri- 
bution of interacting classical particles constrained to the 
presence of an electron on the initial cluster site. Such 
constraint clearly introduces spatial correlations between 
the hopping particle and the environment electrons. By 
creating a "correlation hole" around each carrier, inter- 
actions make polaron hopping in a finite density liquid 
more unfavorable than for non-interacting polarons, im- 
plying a reduction of the mobility. As will be shown in 
the next Section, such static correlations are reflected in 
an increase of the activation barrier for electrical trans- 
port. On the other hand, having implicitly assumed that 
the environment of any given particle is at equilibrium 
(i.e. that it relaxes to equilibrium before the same par- 
ticle can hop again), we are automatically excluding dy- 
namic correlations between subsequent hops. 15 Prelimi- 
nary numerical simulations performed by us on the inter- 
acting liquid indicate that such dynamic correlations can 
at most modify the prefactor of Eq. (I42[) , which amounts 
to logarithmic corrections to the activation barrier. For 
the present problem of interacting polarons, the effect 
would therefore be negligible compared to the effect of 
spatial correlations that we are actually calculating. 

Finally, the textbook result for the mobility of in- 
dependent polarons is recovered by letting £ = in the 
above equations: 



Hp =p 



k B T 



-A P /k B T 



with the prefactor 



UJQ 

t 2 



"I 1/2 



4TA/ 



adiabatic 



non — adiabatic. 



(44) 

(45) 
(46) 
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D. Mean- field approximation 

A complete determination of the statistical distribu- 
tion P(£) defined in the preceding Section requires the 
knowledge of all the many-particle correlation functions 
of the system (generally speaking, the n-th moment of the 
distribution is related to an n-particle correlation func- 
tion). To obtain a tractable expression for the mobility, 
here we evaluate the effect of electron-electron interac- 
tions on the average hopping rate Eq. (|43|) at mean-field 
level, i.e. neglecting the fluctuations of the local field £. 
This scheme of approximation corresponds to the theory 
applied in Ref.^ to the study of organic/dielectric inter- 
faces. It amounts to substituting the averaged hopping 
rate Eq. (|4"3"|) with its first cumulant 



P exp[-A((0)/k B T}. 



(47) 



With this replacement, the problem can be solved in 
terms of the sole two-particle correlation function of the 
interacting system, through the evaluation of the average 
local field 



(c) 



(0 = J2fr2, j -V 1 j]{r 



(48) 



where the symbol (nj)i stands for the constrained prob- 
ability of occupation of site j with site I occupied. As 
anticipated earlier, while the unconstrained average of £ 
would clearly vanish by symmetry in a homogeneous sys- 
tem, the spatial correlations enforced by this constraint 
cause a net additional energy cost (£) > for hopping 
from site to site in the presence of repulsive interactions. 
Assuming that the interaction correction (£) < Ap, so 
that the quadratic term £ 2 in Eq. (|38|) can be neglected, 
we obtain a barrier 



A«0) 



(0 



(49) 



which is the sum of the polaronic activation energy and 
a many-body correction term due to interactions. From 
Eq. (|47p the density dependent mobility can finally be 
expressed in terms of the mobility of independent po- 
larons Eq. ((31)) as 



fi = npexp[-(0/2k B T} 



(50) 



This result shows that in the regime (£) < Ap the many- 
body effects on the mobility arc completely decoupled 
from the individual polaron properties. 

It can be noted that Eqs. (149]) and (|50"|) are formally 
equivalent to the formulas commonly used to describe im- 
purity conduction in compensated polar semiconductors 
and in transition metal oxide glasses i 10 ' 21 ! 24 ' 25 In such 
disordered systems, however, the microscopic mechanism 
responsible for the increase of the polaronic barrier is ex- 
trinsic to the polaronic system, as it originates from the 
ability of the particles to find an efficient percolating path 



connecting dilute, randomly distributed, impurities. 26 
That picture is fundamentally different from the one con- 
sidered here, where (£) originates from the mutual inter- 
actions between carriers in a perfectly crystalline mate- 
rial. 



IV. LONG-RANGE COULOMB INTERACTIONS 

We now apply the theory developed so far to the calcu- 
lation of the mobility of a liquid of small polarons in the 
presence of Coulomb interactions. We shall treat sepa- 
rately the cases of interacting polarons in two and three 
space dimensions: the former applies to the problem of 
polar interfaces as can be found in OFETs with highly 
polarizable gate dielectrics, while the latter can be rel- 
evant for doped polar semiconductors and oxides with 
strong electron-phonon interactions. In both situations, 
the hopping motion associated to the polaronic nature 
of the charge carriers prevents a proper screening of the 
interactions, so that the full long-ranged Coulomb poten- 
tial needs to be considered. We shall therefore take the 
general form 



*\2 



(£) 

Ri 



(51) 



where the effective charge e* accounts for the dielectric 
screening of the polar medium. It is shown in Appendix 
|B1 that e* = ey/2/ (k + e s ) at a two-dimensional polar 
interface, and e* = e/^/eJ in a bulk polar material. 

We start from the observation that for a Coulomb sys- 
tem the correlation function (rij)i appearing in Eq. (|48|) 
varies on lengthscales set by the average inter-particle 
distance ~ R s , defined as R s = (nn)^ 1 ^ 2 in two dimen- 
sions and R s = (47rn/3) -1 / 3 in three dimensions, n being 
the particle density. At sufficiently low concentrations, 
R s is much larger than the lattice spacing so that this 
function can be safely replaced by its continuous limit. 
Correspondingly, the discrete sum appearing in Eq. (|48|) 
can be replaced by the following integral 



(0=n I dr [U(r + R 12 )-U(r)]. 9 ( 2 )(r), 



(52) 



with R\2 — a and g^ 2 \r) the pair distribution function 
of a classical liquid of interacting charged particles - 
the one component plasma (OCP)j2! The properties of 
the OCP are governed by a single dimensionless coupling 
parameter 



r = 



(e*) 2 /Rs 



(53) 



measuring the ratio between the electrostatic interac- 
tions and the thermal energy. This parameter identifies 
a weakly correlated and a strongly correlated regime re- 
spectively for r <C I and r > 1. Upon expanding the 
term between brackets in Eq. (|52p to second order in 
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FIG. 2: The function F(T) for a classical two-dimensional 
Coulomb liquid: Monte-Carlo data of Refl3(J (open circles), 
interpolating function Eq. (|58|l (full line) and Wigner crystal 
result (dashed line). 



y = a/R s , it is readily shown that (£) can be expressed 
in terms of the dimensionless quantities T and y as 



(0 



y 2 F(n 



(54) 



with F(T) a universal function of the OCP. The many- 
body effects on the activation barrier are therefore en- 
tirely controlled by the parameter T characterizing the 
interacting liquid. 



2D 



For a homogeneous two-dimensional system, perform- 
ing the angular integration in Eq. (|52[) an d integrating 
the resulting expression by parts we obtain 



f°° q (2) (v) 

Jo y 



(55) 



In the low density /weakly interacting regime r < 1, 
the correlations of the classical OCP are fully determined 
by the Debye-Hiickel form™ 



(56) 



Upon substituting this function into Eq. (|55[) one obtains 
F{T) = 1. 

In the opposite limit of strong coupling, the electronic 
system undergoes Wigner crystallization, which occurs 
for r > 125j2£ In this regime, it is easy to calculate the 
energy corresponding to a spatial displacement u of a 
given electron while the remaining particles are kept at 
rest. Since the electron under study is initially in an 
equilibrium position, the energy variation is quadratic in 
the displacement and can be written as 



E{u) - E(0) 



; 2m ' 



(57) 



Substituting u — a and converting into the proper units 
we obtain F(T) = (T. The value of the numerical con- 
stant C = 0.8 has been obtained through direct Ewald 
summation of the Coulomb interactions on a triangular 
lattice^ which is the lowest energy structure of a Wigner 
crystal in two dimensions. 

For the evaluation of (£) at intermediate interaction 
strengths we resort to the Monte-Carlo simulations of 
the classical two-dimensional OCP performed in Refl30l. 
There the pair distribution function (r) was tabulated 
at different values of the Coulomb interaction parameter. 
Upon performing the integral Eq. (|55[) using such nu- 
merical data, one obtains a discrete set of points for the 
function F(T). In the range 1 < T < 20, the result can be 
parametrized through the linear interpolating function 



F(T) = 1 + 0.85r 



(58) 



within 1% accuracy (cf. Fig. (2), and this formula re- 
mains fairly accurate even at larger values of T, until it 
eventually merges into the strong coupling Wigner crys- 
tal estimate. It can be observed that, except for a con- 
stant preasymptotic term of order 1, the function F(T) 
representing the interparticle correlations in the Wigner 
crystal has essentially the same T dependence as that of 
the correlated liquid. Using Eqs. |53|) . ([58]) and the defi- 
nition of R s , we can finally write the many-body correc- 
tion to the polaronic activation barrier due to Coulomb 
interactions as 



(0 



fc s r + 0.85(e*) 2 ( 7 rn) 1 / 2 



(59) 



The average local field becomes temperature independent 
and behaves asymptotically as (£) oc n 3 / 2 in the strongly 
correlated limit (r ^> 1), i.e. when the second term be- 
tween brackets dominates. 

Considering the effective polaron-polaron interac- 
tion derived in Appendix [B] for organic/dielectric in- 
terfaces, and using the parameters appropriate to a 
rubrene/Ta 2 5 devices of Ref. 9 (a = 7. 2 A, e s = 25, 
k = 3), we infer that a moderately correlated polaron 
liquid is realized in the two-dimensional conducting chan- 
nel, with coupling parameters in the range < T < 9. 
In this regime polaron-polaron correlations yield an in- 
crease of the activation barrier for transport reaching 
(£)/2 ~ hmeV at the highest concentrations measured. 
This is smaller than the barrier Ap — 55meU for inde- 
pendent polarons, consistent with the assumptions un- 
derlying our derivation. When substituted into Eq. (|50[) . 
such many-body correction leads to a sizable reduction 
of the polaronic mobility, as illustrated in Fig. [3] at two 
different temperatures^ 



B. 3D 

In three space dimensions, 

/>00 7 

F(r) = ry o d y-g^\y)=T. (60) 
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R a = (Trn)" 1 / 2 




fc s T + 0.85(e*) 2 (7rn) 1/2 


3D 


R s = (Ann/3)- 1/:i 


<0 = — (e ) na 



TABLE I: Summary of the main formulas determining the 
density-dependent mobility of small polarons interacting 
through the long-range Coulomb potential of Eq. (|51[) . in two 
and three dimensions. R 3 is the mean interparticle separa- 
tion, a the hopping distance, equal to the distance between 
molecular units, an e* is the effective charge determined by 
the dielectric environment (see Appendix B). The right col- 
umn follows from Eqs. J53J), J54J), J55]) and J60j. It gives the 
mean-field correction to the polaron mobility due to many- 
body effects through fj,/fi P = exp{-{g)/2k B T) [Eg. 1(50)1] ■ 



FIG. 3: Interaction induced reduction of the polaronic 
mobility, calculated with the parameters appropriate to a 
rubrene/Ta2 05 interface, at two different temperatures. The 
full curves are obtained from Eq. (|50[) . while the dashed 
curves include the full activation barrier given by Eq. ((38}, 
with Ap — 55meV . The differences are negligible because in 
all the explored range (£) < Ap. Similar curves are obtained 
for bulk transition-metal oxides. 

This result, which follows directly from the fact that in 
a homogeneous liquid phase <?( 2 )(oo) = 1 (the pair cor- 
relations vanish at large distances), holds exactly at all 
r. It is therefore not necessary to integrate numerically 
the pair distribution function obtained from Monte-Carlo 
simulations as was done in the two-dimensional case. It 
can be directly checked that the result Eq. (|60|) also ex- 
tends to the crystallized phase. To this aim we observe 
that the energy cost to displace a particle from its equi- 
librium position in a three-dimensional Wigner crystal is 
still given by Eq. (f57() . with now ^ 3D ^ = 1 from Gauss' 
theorem^ also leading to F(T) = T. 

For practical calculations the result can be rewritten 
in terms of the carrier density n in a generic three- 
dimensional system as 

(0 = ^-W, (61) 

which is obtained by substituting Eq. ((60)) into Eq. (j5~4")) . 
As usual, a is the hopping distance, of the order of the 
lattice spacing. The linear density dependence of the 
local field resulting from Eq. ((BT]) is weaker than the n 3 / 2 
behavior obtained in two dimensions, and is temperature 
independent at all densities. 

Using typical values for transition-metal oxides such 
as a = 4A and e s = 10 — 100 and assuming a cubic lat- 
tice structure for simplicity we obtain a barrier increase 
(£) /2 = ax, where x is the carrier concentration and the 
coefficient a ~ 40 — AOOmeV. An increase of activation 
energy with electron concentration compatible with such 
prediction has been observed in doped three-dimensional 
transition metal oxides exhibiting small polaron conduc- 
tion, such as magnetite 3 ^ and the manganite a 33 ' 3 ^ . Ac- 
tually, in the high temperature phases of the manganite 



compounds La^Cai-^MnOa, both a large polaron sce- 
nario (in bulk samples^) and a small polaron scenario 
(in thin films^ 3 -) have have been invoked to interpret the 
transport properties in the lightly electron-doped regime. 
We have performed a linear fit of the doping dependence 
of the activation energy A reported in Ref. (33| in the 
range < x < 0.35, yielding A = 46 + 56a; meV. 
When compared with Eqs. ((4*9")) and f6"T)) . the fitted slope 
of the concentration-dependent term yields e s ~ 70, in 
good agreement with the dielectric constants measured 
in those compounds (e s ~ 55 — 90 from Ref. l36l ). 



V. DISCUSSION AND CONCLUSIONS 

In this work we have derived a theory for the hop- 
ping transport of mutually interacting polarons in narrow 
band materials. Observing that in the hopping regime 
the quantum coherence of the carriers extends over only 
few lattice sites, we solve for the quantum dynamics of 
the carriers within a finite size cluster, taking into ac- 
count the interactions with the other charges in the en- 
vironment via a set of static fields. The calculation then 
proceeds by assuming that transport occurs through sta- 
tistically independent hopping events. Correspondingly, 
the many-particle mobility is obtained from a statisti- 
cal average of the inter-molecular hopping rates over the 
distribution of environment fields, which follows from the 
known statistical properties of the interacting liquid. 

The proposed decoupling scheme, which is analogous 
to the one followed by Pardee and Maha n 14 i 15 in the con- 
text of ionic conductors, is justified here by the quasi- 
static nature of the carriers in the hopping regime as a 
consequence of polaronic self-trapping. Despite this sim- 
plification, which amounts to neglecting dynamical corre- 
lations between subsequent hops, the spatial correlations 
between particles which constitute the dominant many- 
body effects on polaronic transport are fully retained. 
When applied to a liquid of small polarons interacting 
through long-range Coulomb forces, the theory predicts 
a net increase of the activation barrier for electrical trans- 
port and hence a reduction of the carrier mobility. The 
analytical formulas obtained at mean-field level, i.e. ne- 
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glecting the fluctuations of the environment field £ repre- 
senting the polaron-polaron correlations, are summarized 
in Table [J 

The present scenario consistently explains the current 
characteristics of rubrene/Ta2C>5 OFETs measured in 
Ref. 1. There, a saturation of the usual linear I oc V g re- 
lationship expected for independent carriers was observed 
at large values of the gate voltage V g , indicative of a siz- 
able reduction of the mobility (the reader is referred to 
that work for a detailed comparison with the experimen- 
tal data). An interpretation in terms of carrier-carrier 
interactions comes naturally in these devices where, as 
was mentioned in the introduction, all the conditions for 
the observation of the predicted many-body effects on 
the polaronic hopping transport are simultaneously met: 
small polaron formation (because of the strong polar cou- 
pling with the gate dielectric and the narrow bandwidth 
of the organic semiconductor), long-range Coulomb re- 
pulsion between the carriers and broad tunability of the 
carrier concentration via the applied gate potential. 

We anticipate based on our theoretical results that, in 
principle, nothing prevents the observation of a down- 
turn of the I — V g curves beyond the saturation regime 
observed in Ref. U For this, the only requirement is 
that of a stronger reduction of the mobility than the one 
realized at rubrene/Ta20s interfaces. As is clear from 
Fig. [21 this can be achieved either by reducing the tem- 
perature, or by increasing the carrier density, as both 
effects lead to an increase of the correlation parameter 
T [see Eq. f55|] and therefore of the ratio (£,)/2k B T in 
Eq. ([50)1 . An interesting possibility in this direction is 
offered by the use of polar electrolytes as gate materials, 
allowing to reach much higher concentrations than with 
conventional polar dielectrics i 37 ' 38 

Finally, due to the very general nature of the mech- 
anisms involved, one might ask if similar effects can be 
observed in other classes of systems. In principle, any 
system with a sufficient concentration of small polarons 
(whatever the microscopic origin) interacting through 
long-range repulsive forces should exhibit a density- 
dependent increase of the transport activation energy. In 
fact, we have found at least two examples in the literature 
which could fit in the present scenario. In the mangan- 
ite compound La^Cai-^MnOa, systematic experimen- 
tal studies of polaronic transport in both thin films 33 
and bulk samples 3 ^ have reported a monotonic increase 
of activation barrier upon increasing the electron con- 
centration x, that could be ascribed to polaron-polaron 
interactions^ A similar increase has been observed in 
Ti doped magnetite (Fc3_j,Tij,C>4) 32 , where a possible 
explanation in terms of long-range Coulomb interactions 
between the carriers has also been explicitly suggested. 
In both classes of compounds, the linear increase of the 
activation energy with electron doping is indeed com- 
patible with the predictions of our theory. Nevertheless, 
other mechanisms can not be excluded, related to the 
complex structural details of these materials, as well as 
to the presence of randomly distributed ionized dopants, 



whose electric fields could also affect the polaronic hop- 
ping rates. 

We conclude by suggesting an experimental method 
that could be useful to disentangle more clearly the ef- 
fects of polaron-polaron interactions from the intrinsic 
features of non- interacting polarons. Such method re- 
lies on the comparison of the activation energy A deter- 
mined from electrical transport, and As obtained from 
thermoelectric power measurements. Since the ther- 
mopower is insensitive to the polaronic renormalization 
of the carriers, As would give a direct measure of the 
interaction correction (£) /2 alone, while electrical trans- 
port would be governed by the sum A = Ap + (£)/2. 
Such method has been often applied to disentangle impu- 
rity effects from polaron effects in transition metal oxide 
glasse s 25 ! 39 and has also been proposed in the context of 
ionic conductors i 15 ! 40 Comparative analysis of the elec- 
trical and thermal transport have also been performed in 
the manganite compounds, to ascertain the polaronic na- 
ture of the charge carriersj 34 i 41 The feasibility of thermo- 
electric power measurements in OFETs has been recently 
demonstrated in Refl42l. and could provide further inde- 
pendent insught into the many-body physics of organic 
field-effect transistors. 

An extension of the present theory to include the ef- 
fects of polaron-polaron correlations beyond the mean- 
field approximation, as well as its generalization to dis- 
ordered systems, is underway. 
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APPENDIX A: HAMILTONIAN OF A TWO SITE 
CLUSTER 

For a two site cluster, the Hamiltonian Eq. ex- 
plicitly reads 

H cl = H t + V h 2n 1 n 2 + Y 1 T] 1 +Y2T] 2 + (Al) 
+ ni(ei + Yx) + n 2 (e 2 + Y 2 ) + 

For the present problem, we can assume without loss of 
generality that the cluster is singly occupied (there is one 
electron on the initial site, the other site being empty for 
the hopping process to be allowed). It is then possible 
to rewrite the cluster Hamiltonian in a form which is 
formally equivalent to that of a tunneling charge inter- 
acting with a single effective mode, which is essentially a 
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spin-boson model. Defining the couplings 
f = [.9 2 ]n + [S 2 ]i2 



g 2 = [.9 2 ]n 



introducing the new variables 



X 



E 



W 2 



Y 2 



(A2) 
(A3) 



Y 2 - Y i 



2fj 

("2 



N = 



V2 
£2 ~ ei 

m + m 
m - m 



V2 



(A4) 



(A5) 



(A6) 



and enforcing the single occupancy within the cluster 



through the condition n\ 
(|AT|) as 



Hi 



E 

71 



V2 



V2 



n 2 = 



1 we can rewrite Eq. 



ni - n 2 
~V2~ 



-kX 2 



-v) + 

+ \kx 2 



(A7) 



From Eq. (|A7|) we see that the variables r\ only contribute 
to an unimportant shift in the x equilibrium position. It 
is therefore convenient to introduce the deviation Q = 
x — grj/k as well as a new interaction variable 



(A8) 



which takes into account both the elecrton-electron in- 
teraction and the electron-phonon screening correction. 

Dropping all terms which do not couple to the site oc- 
cupations or to the phonon displacement, and introduc- 
ing the pseudo-spin notation a z = ri\ — n 2 , H t = —ta x , 
the Hamiltonian H c i can finally be written as 

H sb = ~kQ 2 + -j=(gQ + Z/V2)<T z -ta x . (A9) 



APPENDIX B: EFFECTIVE 
ELECTRON-ELECTRON INTERACTIONS 

1. Organic/dielectric interfaces 

In organic field-effect transistors, charge carriers ac- 
cumulate in a two-dimensional layer located at the in- 
terface between an organic crystal and a polar gate 
dielectric. 1 The model Eq. (Ij therefore consists of two- 
dimensional tight binding electrons interacting with the 
polar phonon modes of the interface. In Fourier space, 



the electron-phonon interaction matrix element has the 
simple forrr>2ii£i^ 



M q = Moe-« z /^q 



(Bl) 



where q is the momentum parallel to the interface, z is 
the distance of the electrons to the polar interface, which 
acts as a short-distance cutoff, and Mq is a coupling con- 
stant that depends on the dielectric properties of the 
interfaced It is given by M§ = 2-Khuioe 2 (3 / 'S , with S 
the total surface of the system, and ujq the frequency of 
the coupled dispersionless polar mode. The parameter 
(3 is a combination of the known dielectric constants of 
the two media that constitute the interface, which de- 
termines the strength of the electron-phonon coupling. 
In the present example of an organic/dielectric interface, 
(3 = (e s — £oo)/(es + k)/(£oo + k) where K is the (fre- 
quency independent) dielectric constant of the organic 
semiconductor, and e s , £ m are respectively the static and 
high-frequency dielectric constants of the polarizable di- 
electric. 

We start with the "bare" interaction potential Vtj be- 
tween two charges located at a distance z from the inter- 
face: 



Vij = - 

K 



1 

Ri, 



1 



Iff 



4z 



2 e c 



(B2) 



where eoo accounts for the high frequency electronic po- 
larizability of the polar material. To determine the effec- 
tive potential we evaluate 



and with Eq. (|2"Tj) we obtain 



e 

K 



1 

Ri: 



■Az 2 ' 



Az 2 
(B3) 

(B4) 



This result is equivalent to what one would obtain from 
a simple image charge calculation, considering the full 
static dielectric constant e s of the polar material right 
from the beginning!^ 

It was shown in Ref. that the conduction in organic 
FETs effectively takes place within the first molecular 
layer nearby the interface. The cut-off distance z is there- 
fore of the order of the lateral size of the molecules, which 
is comparable with the lattice spacing a itself. At concen- 
trations such that the typical inter-particle spacing R s is 
much larger than both a and z, the effective interaction 
potential Eq. (IB4p reduces to 



Vij = 



K Ri 



(B5) 



which corresponds to a long ranged Coulomb potential 
with a screened charge e* = e^Jlj (e s + n). 
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2. Bulk polar materials 



corresponding to a screened charge e* = e/y^ej. 



In three-dimensional polar systems one starts with the 
bare interaction potential 



(B6) 



Local interactions 



where accounts for the high frequency polarizability 
of the material. The interaction of the electrons with 
the polar phonon modes is described by the Frohlich ma- 
trix element M q = M Q /q, with M$ = 27rftw (e 2 /e)/£l 
Here f2 is the total volume of the system, ojq the fre- 
quency of the coupled dispersionless phonon mode and 
e = (e^ 1 — e^ 1 ) -1 an effective dielectric constant. In- 
cluding the screening effect of the polar modes as given 
by Eq. (fT6|) correctly yields 



3 e s R 



(B7) 



To conclude this Appendix we observe that local 
electron-phonon interactions as the ones described by the 
Holstein model do not give rise to a long-range screen- 
ing term. This can be readily seen from Eq. (|2ip . 



where [g 2 ]ij/k oc <Jy. The effective electron-electron in- 
teractions are therefore of the unscreened form Vij = 
e 2 /e QO R i j in bulk materials, and Vij = 2e 2 /(n + e^Rij 
at interfaces. For a given carrier density, the coupling 
parameter T is therefore larger than in the case of po- 
lar screening, and the interaction effects on the mobility 
should be correspondingly enhanced. 
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For the Frohlich interaction in three dimensions a straight- 
forward calculation using the Fourier transform of the ma- 
trix element M q oc 1/q on a cubic lattice gives 7 ~ 0.3, 
while for the electron-phonon interaction at polar inter- 
faces 7 ~ 1/2 — z/y/a 2 + 4z 2 monotonically decreases with 
the distance z to the interface. 

It can also be verified that the condition for adiabatic 
transport is fulfilled in such devices, as the ratio between 
the l.h.s. and the r.h.s. of Eq. (|40|1 is ~ 0.3 at room tem- 
perature. Actually such value places these devices close to 
the adiabatic/nonadiabatic crossover. In this regime, the 
polaronic activation energy changes from Ap = (g 2 /4k) — t 



to the non-adiabatic value Ap = (g 2 /4k), which could be 
at the origin of the missingcorrection —t in the activation 
energy reported in Refs. @S 

It should be stressed that the interaction Eq. (|B1|> was 
derived from the macroscopic laws of electrostatics, that 
are valid at distances > a. In real interfaces, the discrete 
nature of the polarizable medium should lead to an ad- 
ditional short-range cut-off at lengths of the order of the 
inter-ionic spacing. To a first approximation, this effect can 
be incorporated by treating z as an effective phenomeno- 
logical quantity which includes both the channel-interface 
distance and the lattice cutoff. 



